Mathematical Modeling to Explore Shorter Multi-Drug Therapy Options for Active Pulmonary Tuberculosis

Charts for publication

This document details the plots intended for publication of the manuscript “Mathematical Modeling to Explore Shorter Multi-Drug Therapy Options for Active Pulmonary Tuberculosis” primarily authored by John Fors.

Data files were acquired by J Fors from the Systems Pharmacology Model established and implemented in C++ by J Fors.

Figure notation follow that of the final manuscript. All plots were saved as pdf files and imported in Adobe Illustrator for final cleanup and improved legends and exported as .png files at x4 resolution.

Summary of figures

Charts made using this document in sparkline format:

Error in arrangeGrob(...) : object 'Plot1' not found

Libraries & themes

Libraries

library(dplyr)        #data wrangling
library(ggplot2)      #data plotting
library(ggthemes)     #improve plot aesthetics
library(tidyr)        #gather function for multiple columns
library(zoo)          #useful merge and nalocf function
library(gridExtra)    #arrange multiple grid-based plots on a page

• Links to package documents:

dplyr
ggplot2
ggthemes
tidyr
zoo
gridExtra

Porterfield color theme for charts

Color blind friendly theme using dark grey, red, blue, yellow and dark blue

porterfield <- c('#ca222c',  # red
                 '#22a6d6',  # blue
                 '#484848',  # dark grey
                 '#ffca63',  # yellow
                 '#2c7bb6')  # dark blue

Gradient color theme for charts

Blue color gradient for continuous data

colorGradient <- c("#99c5df", #light blue
                   "#3d89c1", #blue
                   "#0a4595", #dark blue
                   "black")   #black  

colorGradient2 <- c("#faae61", #light blue
                   "#ca222c", #blue
                   "#818181", #grey
                   "black")   #black

Original Charts Revisited

Chart 1 - No Treatment

Importing and cleaning data

Data outputs for publication were saved as .csv files and consist of single columns for intracellular (B1) and extracellular (B2) and their respective q1 and q3 when no treatment is on board. Each row is equal to 1 hour.

Reading files:

file1 <- data.frame(B1 = scan('~/Dropbox/John_manuscript/Figure 2/Bt_I_mN.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

file2 <- data.frame(B1q1 = scan('~/Dropbox/John_manuscript/Figure 2/Bt_I_q1N.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

file3 <- data.frame(B1q3 = scan('~/Dropbox/John_manuscript/Figure 2/Bt_I_q3N.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

file4 <- data.frame(B2 = scan('~/Dropbox/John_manuscript/Figure 2/Bt_II_mN.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

file5 <- data.frame(B2q1 = scan('~/Dropbox/John_manuscript/Figure 2/Bt_II_q1N.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

file6 <- data.frame(B2q3 = scan('~/Dropbox/John_manuscript/Figure 2/Bt_II_q3N.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

Merging columns into single data frame and adding time column

a <- bind_cols(file1, file2, file3, file4, file5, file6)
a <- a %>% 
  mutate(time = seq(0, 600, by = 1)) %>% 
  select(time, B1, B1q1, B1q3, B2, B2q1, B2q3)

Final data frame used

a

Final Plot

Graphing the data using ggplot2

Plot1 <- ggplot(a, aes(time, B1)) +
  geom_ribbon(aes(ymin = B1q1, ymax = B1q3), alpha = 0.2, fill = "#21a6d5", size = 0.01) +
  geom_line(color = "#21a6d5") +
    geom_ribbon(aes(ymin = B2q1, ymax = B2q3), alpha = 0.2, fill = "#ca222c", size = 0.01) +
  geom_line(data = a, aes(time, B2), color = "#ca222c") +
  scale_y_log10(breaks = c(1, 10^2, 1e+4, 1e+6, 1e+8)) +
  ylab("Bacterial Load (CFU/mL)") +
  xlab("Time (days since start of infection)") +
  scale_x_continuous(limits = c(0, 600), breaks = seq(0, 600, by = 60)) +
  theme_bw(base_family = "Helvetica") +
  theme(legend.position = 'none')

plot(Plot1)

Figure 2. Bacterial dynamics of simulated acute TB infection. Intra- (red) and extracellular (blue) bacteria count in lung tissue for acute TB infection without drug therapy, for 5000 patients, with infection starting on day 0. Respective shaded areas represent first and third quantiles of population. A bacterial inoculate of 100 CFU/mL infects extracellularly, triggering rapid intra- and extracellular bacterial growth. The acute infection approaches a steady-state level after 9 (270 days) to 12 months (360 days).

Chart 2 - Individual patient

Data outputs for publication were saved as .csv files and data frame consists of 10 columns representing 10 individual patients. Each row is equal to 1 hour.

Reading patient file:

file7 <- read.table('~/Dropbox/John_manuscript/Figure 3/CalcCII2.csv', header = TRUE, sep = ',')

Adding time column and gathering 7 patient columns into single patient column with final wrangling:

file7 <- mutate(file7, time = seq(1, 720, by = 1), days = time/24)
b <- gather(file7, patient, Conc, 1:10)
c <- b %>% 
  group_by(patient) %>% 
  summarise(cmax = max(Conc)) %>% 
  inner_join(b, by ='patient') %>% 
  arrange(cmax)
d <- b %>% 
  group_by(days) %>% 
  summarise(med = median(Conc))

Final data frame

c

d

Final Plot

Graphing the data using ggplot2

Plot2 <- ggplot(c, aes(days, Conc, group = patient, color = patient)) +
  geom_line(alpha = 0.5, size = 0.5, color = "#98AFC7") +
  geom_line(data = d, aes(x = days, y = med, group = NULL),
            colour = "#ca222c", size = 1.5) +
  #scale_y_log10() +                                                    #option to log-scale
  ylab("Concentration (mg/L)") +
  xlab("Time (days since start of treatment)") +
  scale_x_continuous(limits = c(1, 7), breaks = c(1, 2, 3, 4, 5, 6, 7)) +
  theme_bw(base_family = "Helvetica") +
  theme(legend.position = 'none')

plot(Plot2)

Figure 3. Individual patient drug concentration profiles show large variations. Intracellular concentration of rifampin for ten random patients during initial week of drug therapy. Grey lines represent individual patients while the red line represents the median for this group. Large variance is observed, for example, maximum concentration varies more than factor 2x among these patients.

Chart 3 - Drug profiles

Importing and cleaning data

Data outputs for publication were saved as .csv files and data frames consists of single columns representing 4 drugs and their respective intracellular and extracellular data. Each row is equal to 1 hour.

Reading files:

file8 <- data.frame(Conc = scan('~/Dropbox/John_manuscript/Figure 4/concIDrug0_m.txt', skip=1, what = numeric(), sep = '\t'))

Read 721 items

file9 <- data.frame(Conc = scan('~/Dropbox/John_manuscript/Figure 4/concIDrug1_m.txt', skip=1, what = numeric(), sep = '\t'))

Read 721 items

file10 <- data.frame(Conc = scan('~/Dropbox/John_manuscript/Figure 4/concIDrug2_m.txt', skip=1, what = numeric(), sep = '\t'))

Read 721 items

file11 <- data.frame(Conc = scan('~/Dropbox/John_manuscript/Figure 4/concIDrug3_m.txt', skip=1, what = numeric(), sep = '\t'))

Read 721 items

file12 <- data.frame(Conc = scan('~/Dropbox/John_manuscript/Figure 4/concIIDrug0_m.txt', skip=1, what = numeric(), sep = '\t'))

Read 721 items

file13 <- data.frame(Conc = scan('~/Dropbox/John_manuscript/Figure 4/concIIDrug1_m.txt', skip=1, what = numeric(), sep = '\t'))

Read 721 items

file14 <- data.frame(Conc = scan('~/Dropbox/John_manuscript/Figure 4/concIIDrug2_m.txt', skip=1, what = numeric(), sep = '\t'))

Read 721 items

file15 <- data.frame(Conc = scan('~/Dropbox/John_manuscript/Figure 4/concIIDrug3_m.txt', skip=1, what = numeric(), sep = '\t'))

Read 721 items

Adding time column and labeling dataframe drug and cpt.

file8 <- file8 %>% 
  mutate(Drug = "PZA", 
         CPT = "Extra",
         time = seq(0, 720, by = 1), 
         days = time/24)
file9 <- file9 %>% 
  mutate(Drug = "INH", 
         CPT = "Extra",
         time = seq(0, 720, by = 1), 
         days = time/24)
file10 <- file10 %>% 
  mutate(Drug = "RMP", 
         CPT = "Extra",
         time = seq(0, 720, by = 1), 
         days = time/24)
file11 <- file11 %>% 
  mutate(Drug = "EMB", 
         CPT = "Extra",
         time = seq(0, 720, by = 1), 
         days = time/24)
file12 <- file12 %>% 
  mutate(Drug = "PZA", 
         CPT = "Intra",
         time = seq(0, 720, by = 1), 
         days = time/24)
file13 <- file13 %>% 
  mutate(Drug = "INH", 
         CPT = "Intra",
         time = seq(0, 720, by = 1), 
         days = time/24)
file14 <- file14 %>% 
  mutate(Drug = "RMP", 
         CPT = "Intra",
         time = seq(0, 720, by = 1), 
         days = time/24)
file15 <- file15 %>% 
  mutate(Drug = "EMB", 
         CPT = "Intra",
         time = seq(0, 720, by = 1), 
         days = time/24)
e <- bind_rows(file8, file9, file10, file11, file12, file13, file14, file15)

Final data frame

e

Final Plot

Graphing the data using ggplot2

Plot3 <- ggplot(e, aes(days, Conc, color = CPT)) +
  geom_line() +
  #scale_y_log10() +                                                   #option to log-scale
  ylab("Concentraion (mg/L)") +
  xlab("Time (days)") +
  scale_color_manual(values = porterfield) +
  scale_x_continuous(limits = c(0, 14), breaks = c(0, 2,4,6,8,10,12,14)) +
  theme_bw(base_family = "Helvetica") +
  facet_grid(Drug~CPT, scales = "free_y") +
  theme(legend.position = 'none')

plot(Plot3)

Figure 4. Summary drug concentration profiles per drug. Population average concentration of each drug as calculated for extra- (blue) and intracellular (red) compartment, for initial two weeks of drug therapy. This assumes standard therapy and 100% patient adherence.
Note: final graph was recolored and “data junk” removed in Adobe Illustrator

Chart 4 - Standard therapy

Importing and cleaning data

Data outputs for publication were saved as .csv files and consist of single columns for intracellular (B1) and extracellular (B2) and their respective q1 and q3 for treatment on board. Each row is equal to 1 hour.

Reading files:

file22 <- data.frame(AcuteTB = scan('~/Dropbox/John_manuscript/Figure 6/countAcuteTB.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

file23 <- data.frame(ClearedTB = scan('~/Dropbox/John_manuscript/Figure 6/countClearedTB.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

Merging data frames and adding time column.

g <- bind_cols(file22, file23)
g <- mutate(g, time = seq(0, 600, by = 1))

Final data frame

g

Final Plot

Graphing the data using ggplot2

Plot5 <- ggplot(g, aes(time, AcuteTB*100)) +
  geom_line(color = "#21a6d5") +
  geom_line(data = g, aes(time, ClearedTB*100, color = "#ca222c")) +
  #scale_y_log10() +
  ylab("Percentage Population") +
  xlab("Time (days since start of treatment)") +
  scale_x_continuous(limits = c(0, 600), breaks = seq(0, 600, by = 60)) +
  theme_bw(base_family = "Helvetica") +
  theme(legend.position = 'none')

Plot5

Figure 5. Population TB treatment outcome with standard drug therapy. Summary of outcome of median standard drug treatment of acute TB infection, based on simulation of 1000 patients, with standard drug therapy (2HREZ/4HR). The TB infection occurs on day 0, and drug therapy is started on day 180. Acute infection in red color, and cleared status in blue.

Chart 5 - Treatment

Importing and cleaning data

Data outputs for publication were saved as .csv files and consist of single columns for intracellular (B1) and extracellular (B2) and their respective q1 and q3 for treatment on board. Each row is equal to 1 hour.

Reading files:

file16 <- data.frame(B1 = scan("~/Dropbox/John_manuscript/Data/Figure5/Bt_I_m.txt", skip=1, what = numeric(), sep = '\t'))

Read 601 items

file17 <- data.frame(B1q1 = scan('~/Dropbox/John_manuscript/Data/Figure5/Bt_I_q1.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

file18 <- data.frame(B1q3 = scan('~/Dropbox/John_manuscript/Data/Figure5/Bt_I_q3.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

file19 <- data.frame(B2 = scan('~/Dropbox/John_manuscript/Data/Figure5/Bt_II_m.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

file20 <- data.frame(B2q1 = scan('~/Dropbox/John_manuscript/Data/Figure5/Bt_II_q1.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

file21 <- data.frame(B2q3 = scan('~/Dropbox/John_manuscript/Data/Figure5/Bt_II_q3.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

Merging data frames and adding time column.

a2 <- bind_cols(file16, file17, file18, file19, file20, file21)
a2 <- mutate(a2, time = seq(-180, 420, by = 1))

Final data frame

a2

Final Plot

Graphing the data using ggplot2

Plot4 <- ggplot(a2, aes(time, B1)) +
  geom_ribbon(aes(ymin = B1q1, ymax = B1q3), alpha = 0.2, fill = "#21a6d5", size = 0.01) +
  geom_line(color = "#21a6d5") +
  geom_ribbon(aes(ymin = B2q1, ymax = B2q3), alpha = 0.2, fill = "#ca222c", size = 0.01) +
  geom_line(data = a2, aes(time, B2, color = "#ca222c")) +
  scale_y_log10() +
  ylab("Bacterial Load (CFU/mL)") +
  xlab("Time (days since start of treatment)") +
  scale_x_continuous(limits = c(-180, 420), breaks = seq(-180, 420, by = 60)) +
  theme_bw(base_family = "Helvetica") +
  theme(legend.position = 'none')

plot(Plot4)

Figure 6. Intra- vs. extracellular bacterial load in acute TB infection. Intra- and extracellular bacteria count in lung tissue for acute TB infection with standard drug therapy (2HREZ/4HR). The TB infection occurs on day 0, and drug therapy is started on day 180, triggering rapid intra- and extracellular bacterial growth. It is possible bacteria inside infected macrophages form a continued reservoir that may sustain extended infection.

Chart 6 - Incomplete drug therapy

Importing and cleaning data

Data outputs for publication were saved as .csv files and consist of single columns for clearedTB after different treatment conditions. Each row is equal to 1 hour.

Control data imported from previous Figure 6 file. Original file countClearedTB2.txt contains same median as countClearedTB1.txt

Reading files:

file24 <- data.frame(TB1 = scan('~/Dropbox/John_manuscript/Figure 7/countClearedTB1.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

file25 <- data.frame(TB2 = scan('~/Dropbox/John_manuscript/Figure 6/countClearedTB.txt', skip=1, what = numeric(), sep = '\t'))

Read 601 items

Merging data frames and adding time column.

h <- bind_cols(file24, file25)
h <- mutate(h, time = seq(-180, 420, by = 1))

Final data frame

h

Final Plot

Graphing the data using ggplot2

Plot6 <- ggplot(h, aes(time, TB1*100)) +
  geom_line(color = "#21a6d5") +
  geom_line(data = h, aes(time, TB2*100), color = "#ca222c") +
  #scale_y_log10() +                                            #option to log-scale
  ylab("Percentage Population Cured (%)") +
  xlab("Time (days since start of treatment)") +
  scale_x_continuous(limits = c(-180, 420), breaks = seq(-180, 420, by = 60)) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 20)) +
  theme_bw(base_family = "Helvetica") +
  theme(legend.position = 'none')

plot(Plot6)

Figure 7. Reduced successful therapy outcome with incomplete drug therapy. Comparison of patient outcomes for treatment of the median outcome of acute TB infection, with standard drug therapy (2HREZ/4HR) vs. median outcome of therapy excluding isoniazid. Based on simulation of 5000 patients.

Chart 7 - Resistance

Importing and cleaning data

Data outputs for publication were saved as .csv files and consist of single columns for resistant bacteria. Each row is equal to 1 hour.

Reading files:

file26 <- data.frame(Br1_1 = scan('~/Dropbox/John_manuscript/Figure 8/Br1_I_m1.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

file27 <- data.frame(Br1_2 = scan('~/Dropbox/John_manuscript/Figure 8/Br1_II_m1.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

file28 <- data.frame(Br1_3 = scan('~/Dropbox/John_manuscript/Figure 8/Br1_I_m2.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

file29 <- data.frame(Br1_4 = scan('~/Dropbox/John_manuscript/Figure 8/Br1_II_m2.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

Merging data frames and adding time column.

i <- bind_cols(file26, file27, file28, file29)
i <- mutate(i, time = seq(-180, 420, by = 1))

Final data frame

i

Final Plot

Graphing the data using ggplot2

Plot7 <- ggplot(i, aes(time, Br1_2)) +
  geom_line(color = "#21a6d5") +
  geom_line(data = i, aes(time, Br1_4),color = "#ca222c") +
  #scale_y_log10() +
  ylab("Bacterial Load (CFU/ml)") +
  xlab("Time (days since start of treatment)") +
  scale_x_continuous(limits = c(-180, 420), breaks = seq(-180, 420, by = 60)) +
  scale_y_log10() +
  theme_bw(base_family = "Helvetica") +
  theme(legend.position = 'none')

plot(Plot7)

Figure 8. Faster growth of drug-resistant bacterial strains due to incomplete drug therapy. Comparison of growth of rifampin mono-resistant bacterial strain during treatment of acute TB infection, with the median outcome of standard drug therapy (2HREZ/4HR) vs. median outcome of therapy excluding isoniazid. Based on simulation of 5000 patients.

Chart 8 - Drug therapy scenarios

Importing and cleaning data

Data outputs for publication were saved as .csv files and consist of single columns for different treatment outcomes. Each row is equal to 1 hour.

Reading files:

file30 <- data.frame(TB8 = scan('~/Dropbox/John_manuscript/Figure 9/countClearedTB8.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

file31 <- data.frame(TB5 = scan('~/Dropbox/John_manuscript/Figure 9/countClearedTB5.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

file32 <- data.frame(TB4 = scan('~/Dropbox/John_manuscript/Figure 9/countClearedTB4.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

file33 <- data.frame(TB1 = scan('~/Dropbox/John_manuscript/Figure 9/countClearedTB1.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

file34 <- data.frame(TB3 = scan('~/Dropbox/John_manuscript/Figure 9/countClearedTB3.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

Merging data frames and adding time column.

j <- bind_cols(file30, file31, file32, file33, file34)
j <- mutate(j, time = seq(-180, 420, by = 1))

Final data frame

j

Final Plot

Graphing the data using ggplot2

Plot8 <- ggplot(j, aes(time, TB8*100)) +
  geom_line(color = "#21a6d5") +
  geom_line(data = j, aes(time, TB5*100),color = "#ca222c") +
  geom_line(data = j, aes(time + 90, TB4*100),color = "#ffca63") + #time changed to match start of treatment
  geom_line(data = j, aes(time, TB1*100),color = "#484848") +
  geom_line(data = j, aes(time, TB3*100),color = "#2c7bb6") +
  scale_y_continuous(limits = c(0, 100), breaks = c(0, 20, 40, 60, 80, 100)) +
  ylab("Percentage Population Cured (%)") +
  xlab("Time (days since start of treatment)") +
  scale_x_continuous(limits = c(0, 420), breaks = seq(-180, 420, by = 60)) +
  theme_bw(base_family = "Helvetica") +
  theme(legend.position = 'none')

plot(Plot8)

Figure 9. Population TB treatment outcome for different drug therapy scenarios. Comparison of median patient outcomes for treatment of acute TB infection, with standard drug therapy (2HREZ/4HR) vs. therapy with increased rifampin dose (900 mg vs. 600 mg), earlier drug therapy start, shorter drug therapy duration, or less frequent drug dose in continuation phase. Based on simulation of 5000 patients. Note curve with earlier drug therapy start adjusted 90 days for easier comparison.

Note: Figure was cleaned up and colors changed in Adobe Illustrator for final manuscript.

Chart 9 - Patient scenarios

Importing and cleaning data

Data outputs for publication were saved as .csv files and consist of single columns for intracellular (B1) and extracellular (B2) and their respective q1 and q3 for treatment on board. Each row is equal to 1 hour.

Reading files:

file35 <- data.frame(TB7 = scan('~/Dropbox/John_manuscript/Figure 10/countClearedTB7.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

file36 <- data.frame(TB6 = scan('~/Dropbox/John_manuscript/Figure 10/countClearedTB6.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

file37 <- data.frame(TB1 = scan('~/Dropbox/John_manuscript/Figure 10/countClearedTB1.txt', skip = 1, what = numeric(), sep = '\t'))

Read 601 items

Merging data frames and adding time column.

k <- bind_cols(file35, file36, file37)
k <- mutate(k, time = seq(-180, 420, by = 1))

Final data frame

k

Final Plot

Graphing the data using ggplot2

Plot9 <- ggplot(k, aes(time, TB6*100)) +
  geom_line(color = "#21a6d5") +
  geom_line(data = k, aes(time, TB1*100),color = "#ca222c") +
  geom_line(data = k, aes(time, TB7*100),color = "#ffca63") +
  scale_y_continuous(limits = c(0, 100), breaks = c(0, 20, 40, 60, 80, 100)) +
  ylab("Percentage Population Cured (%)") +
  xlab("Time (days since start of treatment)") +
  scale_x_continuous(limits = c(0, 420), breaks = seq(0, 420, by = 60)) +
  theme_bw(base_family = "Helvetica") +
  theme(legend.position = 'none')

plot(Plot9)

Figure 10. Population TB treatment outcome for different patient scenarios. Comparison of median patient outcomes for treatment of acute TB infection, with standard drug therapy (2HREZ/4HR), for patients with reduced adherence level (average 0.880% vs. default 100%), or with reduced immune system function (0.550% of normal). Based on simulation of 5000 patients. Infection occurs on day 0, and drug therapy is started at day 180.

New Added Charts

A sensitivity analysis was performed on the model by William Fox using the C++ model. Variables will scaled by 0.1, 0.5, 2, and 10, and 10 iteration performed. The data for the final median acute TB was extracted and saved as .csv files.

Chart 10 - Drug Tornado plot

Importing and cleaning data

The original .csv file was altered in Excel and removed the first row number column and renamed the CLEARED.TB column to omit the full stop. The column name MULTIPLIER was changed to parameterScaleFactor.

The name CLEAREDTB was an artefact of the original sensitivity anlysis. The data is Acute TB, but the column name was not changed.

Reading files:

file38 <- read.csv("sensitivityA1.csv")

Data wrangling to obtain the median of the 10 iterations and filter data frame to individual drug data frames

tornadoMedians <- file38 %>% 
  group_by(DRUG, VARIABLE, parameterScaleFactor) %>% 
  summarise(MEDIAN = median(CLEAREDTB)) %>% 
  filter(parameterScaleFactor != 1) %>%  
          filter(VARIABLE != "decayFactor") %>%
         filter(VARIABLE != "killIntra") %>%
         filter(VARIABLE != "killExtra") %>%
         filter(VARIABLE != "growIntra") %>%
         filter(VARIABLE != "growExtra") %>%
         filter(VARIABLE != "factorFast") %>%
         filter(VARIABLE != "factorSlow")
tornadoMedians$parameterScaleFactor <- as.factor(tornadoMedians$parameterScaleFactor)
#Rename VARIABLES
tornadoMedians[tornadoMedians$VARIABLE == "KaMean", "NAME"] = "Ka (hr-1)"
tornadoMedians[tornadoMedians$VARIABLE == "KeMean", "NAME"] = "Ke (hr-1)"
tornadoMedians[tornadoMedians$VARIABLE == "V1Mean", "NAME"] = "V/F (L)"
tornadoMedians[tornadoMedians$VARIABLE == "KeMult", "NAME"] = "auto-induction factor"
tornadoMedians[tornadoMedians$VARIABLE == "KeTime", "NAME"] = "Time to max auto-induction (d)"
tornadoMedians[tornadoMedians$VARIABLE == "highAcetFactor", "NAME"] = "Clearance multiplier for fast acetylators"
tornadoMedians[tornadoMedians$VARIABLE == "IOfactor", "NAME"] = "Intracellular macrophage Conc multiplier"
tornadoMedians[tornadoMedians$VARIABLE == "V1Mean", "NAME"] = "V/F (L)"
tornadoMedians[tornadoMedians$VARIABLE == "GRfactor", "NAME"] = "Granuloma concentration multiplier"
tornadoMedians[tornadoMedians$VARIABLE == "V1Mean", "NAME"] = "V/F (L)"
tornadoMedians[tornadoMedians$VARIABLE == "EC50k", "NAME"] = "Bacterial killing EC50 (mg/L)"
tornadoMedians[tornadoMedians$VARIABLE == "EC50g", "NAME"] = "Bacterial growth inhibition EC50 (mg/L)"
tornadoMedians[tornadoMedians$VARIABLE == "ak", "NAME"] = "Bacterial killing Hill factor"
tornadoMedians[tornadoMedians$VARIABLE == "ag", "NAME"] = "Bacterial growth inhibition Hill factor"
tornadoMedians[tornadoMedians$VARIABLE == "ak", "NAME"] = "Bacterial killing Hill factor"
tornadoMedians[tornadoMedians$VARIABLE == "mutationRate", "NAME"] = "Bacterial mutation rate (h-1)"
tornadoMedians[tornadoMedians$VARIABLE == "kill_e0", "NAME"] = "Extracellular kill rate days < 0  (h-1)"
tornadoMedians[tornadoMedians$VARIABLE == "kill_e1", "NAME"] = "Extracellular kill rate days 0-2  (h-1)"
tornadoMedians[tornadoMedians$VARIABLE == "kill_e2", "NAME"] = "Extracellular kill rate days 3-15  (h-1)"
tornadoMedians[tornadoMedians$VARIABLE == "kill_e3", "NAME"] = "Extracellular kill rate days 15+  (h-1)"
tornadoMedians[tornadoMedians$VARIABLE == "kill_i0", "NAME"] = "Intracellular kill rate days < 0  (h-1)"
tornadoMedians[tornadoMedians$VARIABLE == "kill_i1", "NAME"] = "Intracellular kill rate days 0-2  (h-1)"
tornadoMedians[tornadoMedians$VARIABLE == "kill_i2", "NAME"] = "Intracellular kill rate days 3-15  (h-1)"
tornadoMedians[tornadoMedians$VARIABLE == "kill_i3", "NAME"] = "Intracellular kill rate days 15+  (h-1)"
tornadoMedians[tornadoMedians$VARIABLE == "multiplier", "NAME"] = "Granuloma multiplier"
tornadoMedians[tornadoMedians$VARIABLE == "rise50", "NAME"] = "Granuloma time constant in"
tornadoMedians[tornadoMedians$VARIABLE == "fall50", "NAME"] = "Granuloma time constant in"
R <- filter(tornadoMedians, DRUG == "RIF", VARIABLE != "highAcetFactor"
            )
H <- filter(tornadoMedians, DRUG == "INH", VARIABLE != "KeMult")
Z <- filter(tornadoMedians, DRUG == "PZA", VARIABLE != "highAcetFactor", VARIABLE != "KeMult")
E <- filter(tornadoMedians, DRUG == "EMB", VARIABLE != "highAcetFactor", VARIABLE != "KeMult")

RIF data frame

R

H

PZA data frame

Z

EMB data frame

E

Building Final Plot

Function to build consistent plot for all drugs where x = df and y = “df Title.” Plot uses change in percentage points

mytornadoP <- function(x, y){
  ggplot(x, aes(x = reorder(NAME, MEDIAN*100-5), y = MEDIAN*100-5, color = parameterScaleFactor)) +
  geom_bar(position = "identity", stat = "identity", fill = "white") +
  scale_color_manual(values = colorGradient2) +
  coord_flip() +
  ggtitle(y) +
  theme_minimal(base_family = "Helvetica") +
  xlab("Parameter") +
  ylab("Change in Acute TB") +
  theme(legend.position = "none")
}

Function to build consistent plot for all drugs where x = df and y = “df Title.” Plot uses fold change

mytornadoF <- function(x, y){
  ggplot(x, aes(x = reorder(NAME, MEDIAN), y = log2(MEDIAN/0.05), color = parameterScaleFactor)) +
  geom_bar(position = "identity", stat = "identity", fill = "white") +
  scale_color_manual(values = colorGradient2) +
  coord_flip() +
  ggtitle(y) +
  theme_minimal(base_family = "Helvetica") +
  xlab("Parameter") +
  ylab("Median log2 fold change") +
  theme(legend.position = "none")
}

Individual change in percentage points plots

#Plot for RIF
PlotR <- mytornadoP(R, "RIF") +
  scale_y_continuous(limits = c(-5, 100), breaks = c(-5, 0, 20, 40, 60, 80, 100))

Error in mytornadoP(R, "RIF") : could not find function "mytornadoP"

Individual fold change plots

#Plot for RIF
PlotRf <- mytornadoF(R, "RIF") 
#Plot for INH
PlotHf <- mytornadoF(H, "INH") 
#Plot for PZA
PlotZf <- mytornadoF(Z, "PZA")
#Plot for EMB
PlotEf <- mytornadoF(E, "EMB")

Final Plot

Plot arranged using package gridExtra. Additional plot using legend was saved for clean up in Adobe Illustrator

grid.arrange(PlotR, PlotH, PlotZ, PlotE)

Figure TBD. . TBD. Note: Plot was ouput and saved as PDF using quartz() to obatin vector file that was scalable and easier to improve aesthetics in Adobe Illustrator

grid.arrange(PlotRf, PlotHf, PlotZf, PlotEf)

Chart 11 - Init Tornado Plot

Importing and cleaning data

The original .csv file was altered in Excel and removed the first row number column and renamed the CLEARED.TB column to omit the full stop. The column name MULTIPLIER was changed to parameterScaleFactor.

The name CLEAREDTB was an artefact of the original sensitivity anlysis. The data is Acute TB, but the column name was not changed.

Read data:

file39 <- read.csv("sensitivityOther.csv")

Data wrangling to obtain the median of the 10 iterations

tornadoMediansO <- file39 %>% 
  group_by(DRUG, VARIABLE, parameterScaleFactor) %>% 
  summarise(MEDIAN = median(CLEAREDTB)) %>% 
  filter(parameterScaleFactor != 1, 
         VARIABLE != "adherenceStdv",
         VARIABLE != "adherenceStdvDay",
         VARIABLE != "adherenceSwitchDay",
         VARIABLE != "bactThreshold",
         VARIABLE != "bactThresholdRes",
         VARIABLE != "growthLimit",
         VARIABLE != "immuneStdv",
         VARIABLE != "infI",
         VARIABLE != "infII",
         VARIABLE != "infIII",
         VARIABLE != "infIV",
         VARIABLE != "initialValueStdv",
         VARIABLE != "parameterStdv",
         VARIABLE != "persistTime",
         VARIABLE != "nTime",
         VARIABLE != "nPatients",
         VARIABLE != "timeStepStdv",
         VARIABLE != "nIterations",
         VARIABLE != "nPopulations",
         VARIABLE != "shareLowAcetylators",
         VARIABLE != "therapyStart")
tornadoMediansO$parameterScaleFactor <- as.factor(tornadoMediansO$parameterScaleFactor)

tornadoMediansO

write.csv(tornadoMediansO, file = "tornadoMediansO.csv")

Extra data wrangling from John Fors done in Excel. Included 1. added full parameter definitions under column name, 2. removed 10 and 5 scale factor for parameters that were ratios/level [0..1]

tornadoMedians1 <- read.csv("tornadoMedians1.csv")
tornadoMedians1$parameterScaleFactor <- as.factor(tornadoMedians1$parameterScaleFactor)

Final Plot

mytornado2 <- function(x, y){
  ggplot(x, aes(x = reorder(VARIABLE, MEDIAN), y = MEDIAN*100-5, color = parameterScaleFactor)) +
  geom_bar(position = "identity", stat = "identity", fill = "white") +
  scale_color_manual(values = colorGradient2) +
  coord_flip() +
  ggtitle(y) +
  theme_minimal(base_family = "Helvetica") +
  xlab("Parameter") +
  ylab("Change in Acute TB") +
  theme(legend.position = "none")
}

mytornado2f <- function(x, y){
  ggplot(x, aes(x = reorder(Name, MEDIAN), y = log2(MEDIAN/0.05), color = parameterScaleFactor)) +
  geom_bar(position = "identity", stat = "identity", fill = "white") +
  scale_color_manual(values = colorGradient2) +
  coord_flip() +
  ggtitle(y) +
  theme_minimal(base_family = "Helvetica") +
  xlab("Parameter") +
  ylab("Median log2 fold change") +
  theme(legend.position = "none")
}

Plot0 <- mytornado2(tornadoMediansO, "Init parameters") +
  scale_y_continuous(limits = c(-10, 100), breaks = c(-10, 0, 20, 40, 60, 80, 100))

Plot0f <- mytornado2f(tornadoMedians1, "Init parameters")

plot(Plot0)

Figure TBD. . TBD. Note: Plot was ouput and saved as PDF using quartz() to obatin vector file that was scalable and easier to improve aesthetics in Adobe Illustrator

plot(Plot0f)

Chart 12 - Immune Tornado plot

Importing and cleaning data

The original .csv file was altered in Excel and removed the first row number column. The column name MULTIPLIER was changed to parameterScaleFactor. The column DRUG was renamed Immune and used to rename parameters according to their specific immune compartment effect, namely: “Bacteria” “Cytokines” “Dendritic_cells” “DLN_IL12” “DLN_Tcell” “Macrophage” “Tcells” as stated by the setImmuneParameters.cpp file in the Source files of the Systems Model C++ tool. This was done in Excel as both T cells and Bacteria compartments had a parameter named “a2”. In excel each iteration was saved in order and was manually changed in sequence rather than using R which grouped the a2 parameters as a single factor.

Reading files:

file40 <- read.csv("sensitivity_I.csv")

Data wrangling to obtain the median of the 10 iterations and filter data frame to individual drug data frames

tornadoMediansI <- file40 %>% 
  group_by(Immune, VARIABLE, parameterScaleFactor) %>% 
  summarise(MEDIAN = median(MedianAcuteTB)) %>% 
  filter(parameterScaleFactor != 1)
tornadoMedians$parameterScaleFactor <- as.factor(tornadoMedians$parameterScaleFactor)
B <- filter(tornadoMediansI, Immune == "Bacteria")
C <- filter(tornadoMediansI, Immune == "Cytokines")
D <- filter(tornadoMediansI, Immune == "Dendritic_cells")
DT <- filter(tornadoMediansI, Immune == "DLN_Tcell")
DI <- filter(tornadoMediansI, Immune == "DLN_IL12")
M <- filter(tornadoMediansI, Immune == "Macrophage")
Tc <- filter(tornadoMediansI, Immune == "Tcells")
B$parameterScaleFactor <- as.factor(B$parameterScaleFactor)
C$parameterScaleFactor <- as.factor(C$parameterScaleFactor)
D$parameterScaleFactor <- as.factor(D$parameterScaleFactor)
DT$parameterScaleFactor <- as.factor(DT$parameterScaleFactor)
DI$parameterScaleFactor <- as.factor(DI$parameterScaleFactor)
M$parameterScaleFactor <- as.factor(M$parameterScaleFactor)
Tc$parameterScaleFactor <- as.factor(Tc$parameterScaleFactor)
DLN <- bind_rows(DT, DI)

write.csv(B, "B.csv")
write.csv(C, "C.csv")
write.csv(D, "D.csv")
write.csv(DT, "DT.csv")
write.csv(DI, "DI.csv")
write.csv(Tc, "Tc.csv")
write.csv(M, "M.csv")

Files converted to .csv and imported to Excel to add Variable definition (time contraints)

B1 <- read.csv("B1.csv")
C1 <- read.csv("C1.csv")
D1 <- read.csv("D1.csv")
DT1 <- read.csv("DT1.csv")
DI1 <- read.csv("DI1.csv")
Tc1 <- read.csv("Tc1.csv")
M1 <- read.csv("M1.csv")
DLN1 <- read.csv("DLN.csv")
B1$parameterScaleFactor <- as.factor(B1$parameterScaleFactor)
C1$parameterScaleFactor <- as.factor(C1$parameterScaleFactor)
D1$parameterScaleFactor <- as.factor(D1$parameterScaleFactor)
DT1$parameterScaleFactor <- as.factor(DT1$parameterScaleFactor)
DI1$parameterScaleFactor <- as.factor(DI1$parameterScaleFactor)
Tc1$parameterScaleFactor <- as.factor(Tc1$parameterScaleFactor)
M1$parameterScaleFactor <- as.factor(M1$parameterScaleFactor)
DLN1$parameterScaleFactor <- as.factor(DLN1$parameterScaleFactor)

Bacteria data frame

B

Cytokine data frame

C

Dendritic cells data frame

D

DLN data frame

DLN

Macrophage cells data frame

M

T cells data frame

Tc

Building Final Plot

Function:mytornado to build consistent plot was used for immune parameters

mytornadoIm <- function(x, y){
  ggplot(x, aes(x = reorder(VARIABLE, MEDIAN*100-5), y = MEDIAN*100-5, color = parameterScaleFactor)) +
  geom_bar(position = "identity", stat = "identity", fill = "white") +
  scale_color_manual(values = colorGradient2) +
  coord_flip() +
  ggtitle(y) +
  theme_minimal(base_family = "Helvetica") +
  xlab("Parameter") +
  ylab("Change in Acute TB") +
  theme(legend.position = "none")
}

mytornadoImf <- function(x, y){
  ggplot(x, aes(x = reorder(Name, MEDIAN), y = log2(MEDIAN/0.05), color = parameterScaleFactor)) +
  geom_bar(position = "identity", stat = "identity", fill = "white") +
  scale_color_manual(values = colorGradient2) +
  coord_flip() +
  ggtitle(y) +
  theme_minimal(base_family = "Helvetica") +
  xlab("Parameter") +
  ylab("Median log2 fold change") +
  theme(legend.position = "none")
}

Individual plots

#Plot for Bacteria
PlotB <- mytornadoIm(B, "Bacteria") +
  scale_y_continuous(limits = c(-5, 100), breaks = c(-5, 0, 20, 40, 60, 80, 100))
#Plot for Cytokines
PlotC <- mytornadoIm(C, "Cytokines") +
  scale_y_continuous(limits = c(-5, 70), breaks = c(-5, 0, 20, 40, 60))
#Plot for dendritic cell
PlotD <- mytornadoIm(D, "Dendritic cells") +
  scale_y_continuous(limits = c(-5, 25), breaks = c(-5, 0, 5, 10, 15, 20, 25))
#Plot for DLN
PlotDLN <- mytornadoIm(DLN, "DLN") +
  scale_y_continuous(limits = c(-5, 32), breaks = c(-5, 0, 10, 20, 30))
#Plot for macrophage cells
PlotM <- mytornadoIm(M, "Macrophages") +
  scale_y_continuous(limits = c(-5, 68), breaks = c(-5, 0, 20, 40, 60))
#Plot for T cells
PlotTc <- mytornadoIm(Tc, "T cells") +
  scale_y_continuous(limits = c(-5, 85), breaks = c(-5, 0, 20, 40, 60, 80))

#Plot for Bacteria
PlotBf <- mytornadoImf(B1, "Bacteria")
#Plot for Cytokines
PlotCf <- mytornadoImf(C1, "Cytokines")
#Plot for dendritic cell
PlotDf <- mytornadoImf(D1, "Dendritic cells")
#Plot for DLN
PlotDLNf <- mytornadoImf(DLN1, "DLN")
#Plot for macrophage cells
PlotMf <- mytornadoImf(M1, "Macrophages")
#Plot for T cells
PlotTcf <- mytornadoImf(Tc1, "T cells")

Final Plot

Plot arranged using package gridExtra. Additional plot using legend was saved for clean up in Adobe Illustrator

grid.arrange(PlotB, PlotC, PlotD, PlotDLN, PlotM, PlotTc)

Figure TBD. . TBD. Note: Plot was ouput and saved as PDF using quartz() to obtain vector file that was scalable and easier to improve aesthetics in Adobe Illustrator

grid.arrange(PlotBf, PlotCf, PlotDf, PlotDLNf, PlotMf, PlotTcf)

PlotTcf

Error in reorder(Name, MEDIAN) : object 'Name' not found